TPC/Upgrade/macros/spaceChargeFluctuation.C

   1 /*
   2   .x ~/rootlogon.C
   3   .L $ALICE_ROOT/TPC/Upgrade/macros/spaceChargeFluctuation.C+
   4 //  GenerateMapRaw("/hera/alice/local/filtered/alice/data/2010/LHC10e/000129527/raw/merged_highpt_1.root","output.root", 50);
   5   GenerateMapRawIons( "/hera/alice/local/filtered/alice/data/2010/LHC10e/000129527/raw/merged_highpt_1.root","output.root", 25);
   6
   7  */
   8 #include "TMath.h"
   9 #include "TRandom.h"
  10 #include "TTreeStream.h"
  11 #include "TVectorD.h"
  12 #include "TCanvas.h"
  13 #include "TStopwatch.h"
  14 #include "AliTPCParam.h"
  15 #include "AliTPCcalibDB.h"
  16 #include "AliTPCAltroMapping.h"
  17 #include "AliAltroRawStream.h"
  18 #include "AliSysInfo.h"
  19 #include "AliTPCRawStreamV3.h"
  20 #include "AliCDBManager.h"
  21 #include "TGeoGlobalMagField.h"
  22 #include "AliMagF.h"
  23 #include "AliRawReaderRoot.h"
  24 #include "AliRawReader.h"
  25 #include "TH3.h"
  26 #include "TH2.h"
  27
  28
  29 TH1 * GenerateMapRawIons(const char *fileName="raw.root", const char *outputName="histo.root", Int_t maxEvents=25);
  30
  31
  32 void spaceChargeFluctuation(Int_t mode=0){
  33   //
  34   //
  35   if (mode==0) GenerateMapRawIons();
  36
  37 }
  38
  39 void spaceChargeFluctuationMC(Int_t nframes){
  40   //
  41   // Toy MC to generate space charge fluctuation
  42   //
  43   //
  44   TTreeSRedirector *pcstream = new TTreeSRedirector("spaceChargeFluctuation.root","recreate");
  45   Double_t interactionRate=100000;
  46   Double_t driftTime=0.1;
  47
  48   for (Int_t iframe=0; iframe<nframes; iframe++){
  49     Int_t nevents=gRandom->Poisson(interactionRate*driftTime);
  50     Int_t ntracksAll=0;
  51     TVectorD vecAllPhi128(128);
  52     TVectorD vecAllPhi36(36);
  53     for (Int_t ievent=0; ievent<nevents; ievent++){
  54       //
  55       //
  56       Int_t ntracks=gRandom->Exp(500)+gRandom->Exp(500); // to be taken form the distribution
  57       ntracksAll+=ntracks;
  58       for (Int_t isec=0; isec<128; isec++){
  59         vecAllPhi128(isec)+=ntracks*gRandom->Rndm()/128;
  60       }
  61       for (Int_t isec=0; isec<36; isec++){
  62         vecAllPhi36(isec)+=ntracks*gRandom->Rndm()/36;
  63       }
  64     }
  65     (*pcstream)<<"ntracks"<<
  66       "nevents="<<nevents<<                  // number of events withing time frame
  67       "ntracksAll="<<ntracksAll<<            // number of tracks within  time frame
  68       "vecAllPhi36="<<&vecAllPhi36<<         // number of tracks in phi bin (36 bins)    within time frame
  69       "vecAllPhi128="<<&vecAllPhi128<<       // number of tracks in phi bin (128 bins)   within time frame
  70       "\n";
  71   }
  72   delete pcstream;
  73 }
  74
  75
  76 // void spaceChargeFluctuationDraw(){
  77 //   //
  78 //   //
  79 //   //
  80 //   TFile *f = TFile::Open("spaceChargeFluctuation.root");
  81 //   TTree * tree = (TTree*)f->Get("ntracks");
  82 //   TCanvas * canvasFluc = new TCanvas("cavasFluc","canvasFluc");
  83 // }
  84
  85
  86 TH1 * GenerateMapRaw(const char *fileName, const char *outputName, Int_t maxEvents){
  87   //
  88   // Generate 3D maps of electrons
  89   // different threshold considered
  90   // Paramaters:
  91   //    fileName      - name of input raw file
  92   //    outputName    - name of output file with the space charge histograms
  93   //    maxEvents     - grouping of the events
  94   //
  95   //
  96   TTreeSRedirector * pcstream  = new TTreeSRedirector(outputName, "recreate");
  97
  98   const char *  ocdbpath =gSystem->Getenv("OCDB_PATH") ? gSystem->Getenv("OCDB_PATH"):"local://$ALICE_ROOT/OCDB/";
  99   AliCDBManager * man = AliCDBManager::Instance();
 100   man->SetDefaultStorage(ocdbpath);
 101   man->SetRun(0);
 102   TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG,       AliMagF::kBeamTypepp, 2.76/2.));
 103   AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
 104   AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
 105   //
 106   TStopwatch timer;
 107   timer.Start();
 108   TH3F * hisQ3D[3]={0};
 109   TH2F * hisQ2D[3]={0};
 110   TH1F * hisQ1D[3]={0};
 111   //
 112   for (Int_t ith=0; ith<3; ith++){
 113     char chname[100];
 114     snprintf(chname,100,"hisQ3D_Th%d",2*ith+2);
 115     hisQ3D[ith] = new TH3F(chname, chname,180, 0,2*TMath::TwoPi(),param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) ,100,0,1100);
 116     snprintf(chname,100,"hisQ2D_Th%d",2*ith+2);
 117     hisQ2D[ith] = new TH2F(chname,chname,180, 0,2*TMath::TwoPi(), param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
 118     //
 119     snprintf(chname,100,"hisQ1D_Th%d",2*ith+2);
 120     hisQ1D[ith] = new TH1F(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
 121     hisQ3D[ith]->SetDirectory(0);
 122     hisQ1D[ith]->SetDirectory(0);
 123     hisQ1D[ith]->SetDirectory(0);
 124   }
 125   //
 126   //
 127   AliRawReader *reader = new AliRawReaderRoot(fileName);
 128   reader->Reset();
 129   AliAltroRawStream* stream = new AliAltroRawStream(reader);
 130   stream->SelectRawData("TPC");
 131   Int_t evtnr=0;
 132   Int_t chunkNr=0;
 133   //
 134   while (reader->NextEvent()) {
 135     if(evtnr>maxEvents) {
 136       chunkNr++;
 137       evtnr=0;
 138       pcstream->GetFile()->mkdir(Form("Chunk%d",chunkNr));
 139       pcstream->GetFile()->cd(Form("Chunk%d",chunkNr));
 140       for (Int_t ith=0; ith<3; ith++){
 141         hisQ1D[ith]->Write(Form("His1D_%d",ith));
 142         hisQ2D[ith]->Write(Form("His2D_%d",ith));
 143         hisQ3D[ith]->Write(Form("His3D_%d",ith));
 144         if (ith<3) continue;
 145         hisQ1D[ith]->Reset();
 146         hisQ2D[ith]->Reset();
 147         hisQ3D[ith]->Reset();
 148       }
 149     }
 150     cout<<"Chunk=\t"<<chunkNr<<"\tEvt=\t"<<evtnr<<endl;
 151     evtnr++;
 152     AliSysInfo::AddStamp(Form("Event%d",evtnr),evtnr);
 153     AliTPCRawStreamV3 input(reader,(AliAltroMapping**)mapping);
 154     //
 155     while (input.NextDDL()){
 156       Int_t sector = input.GetSector();
 157       while ( input.NextChannel() ) {
 158         Int_t    row    = input.GetRow();
 159         Int_t rowAbs    = row+ ((sector<36) ? 0: param->GetNRow(0));
 160         Int_t    pad    = input.GetPad();
 161         Int_t    nPads   = param->GetNPads(sector,row);
 162         Double_t localX  = param->GetPadRowRadii(sector,row);
 163         Double_t localY  = (pad-nPads)*param->GetPadPitchWidth(sector);
 164         Double_t localPhi= TMath::ATan2(localY,localX);
 165         Double_t phi     = TMath::Pi()*((sector%18)+0.5)/9+localPhi;
 166         if (sector%36>=18) phi+=TMath::TwoPi();
 167         Double_t padLength=param->GetPadPitchLength(sector,row);
 168         //      Double_t padWidth=param->GetPadPitchWidth(sector);
 169         //    Double_t zLength=param->GetZLength();
 170         Double_t deltaPhi=hisQ3D[0]->GetXaxis()->GetBinWidth(0);
 171         Double_t deltaZ= hisQ3D[0]->GetZaxis()->GetBinWidth(0);
 172         //
 173         while ( input.NextBunch() ){
 174           Int_t  startTbin    = (Int_t)input.GetStartTimeBin();
 175           Int_t  bunchlength  = (Int_t)input.GetBunchLength();
 176           const UShort_t *sig = input.GetSignals();
 177           Int_t aboveTh[3]={0};
 178           for (Int_t i=0; i<bunchlength; i++){
 179             for (Int_t ith=0; ith<3; ith++){
 180               if (sig[i]>(ith*2)+2) aboveTh[ith]++;
 181             }
 182           }
 183           for (Int_t ith=0; ith<3; ith++){
 184             if (aboveTh[ith%3]>1){
 185               for (Int_t i=0; i<bunchlength; i++){
 186                 //
 187                 // normalization
 188                 //
 189                 Double_t volume3D = padLength*(localX*deltaPhi)*deltaZ;
 190                 if (sector%36<18) hisQ1D[ith]->Fill(rowAbs, sig[i]/volume3D);
 191                 hisQ2D[ith]->Fill(phi,rowAbs, sig[i]/volume3D);
 192                 hisQ3D[ith]->Fill(phi,rowAbs,startTbin+i,sig[i]/volume3D);
 193               }
 194             }
 195           }
 196         }
 197       }
 198     }
 199   }
 200   timer.Print();
 201   delete pcstream;
 202   return 0;
 203 }
 204
 205
 206
 207
 208
 209 TH1 * GenerateMapRawIons(const char *fileName, const char *outputName, Int_t maxEvents){
 210   //
 211   // Generate 3D maps of the space charge for the rad data maps
 212   // different threshold considered
 213   // Paramaters:
 214   //    fileName      - name of input raw file
 215   //    outputName    - name of output file with the space charge histograms
 216   //    maxEvents     - grouping of the events
 217   //
 218   //
 219   TTreeSRedirector * pcstream  = new TTreeSRedirector(outputName, "recreate");
 220   const char *  ocdbpath =gSystem->Getenv("OCDB_PATH") ? gSystem->Getenv("OCDB_PATH"):"local://$ALICE_ROOT/OCDB/";
 221   AliCDBManager * man = AliCDBManager::Instance();
 222   man->SetDefaultStorage(ocdbpath);
 223   man->SetRun(0);
 224   TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG,       AliMagF::kBeamTypepp, 2.76/2.));
 225   AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
 226   AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
 227   TStopwatch timer;
 228   timer.Start();
 229   //   arrays of space charges - different elements corresponds to different threshold to accumulate charge
 230   TH3F * hisQ3D[3]={0};                // 3D maps space charge from drift volume
 231   TH2F * hisQ2D[3]={0};
 232   TH1F * hisQ1D[3]={0};
 233   TH3F * hisQ3DROC[3]={0};             // 3D maps space charge from ROC
 234   TH2F * hisQ2DROC[3]={0};
 235   TH1F * hisQ1DROC[3]={0};
 236   //
 237   for (Int_t ith=0; ith<3; ith++){
 238     char chname[100];
 239     snprintf(chname,100,"hisQ3D_Th%d",2*ith+2);
 240     hisQ3D[ith] = new TH3F(chname, chname,180, 0,2*TMath::TwoPi(),param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) ,125,-250,250);
 241     snprintf(chname,100,"hisQ2D_Th%d",2*ith+2);
 242     hisQ2D[ith] = new TH2F(chname,chname,180, 0,2*TMath::TwoPi(), param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
 243     snprintf(chname,100,"hisQ1D_Th%d",2*ith+2);
 244     hisQ1D[ith] = new TH1F(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
 245     //
 246     snprintf(chname,100,"hisQ3DROC_Th%d",2*ith+2);
 247     hisQ3DROC[ith] = new TH3F(chname, chname,180, 0,2*TMath::TwoPi(),param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) ,125,-250,250);
 248     snprintf(chname,100,"hisQ2DROC_Th%d",2*ith+2);
 249     hisQ2DROC[ith] = new TH2F(chname,chname,180, 0,2*TMath::TwoPi(), param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
 250     snprintf(chname,100,"hisQ1DROC_Th%d",2*ith+2);
 251     hisQ1DROC[ith] = new TH1F(chname,chname, param->GetNRow(0)+param->GetNRow(36) ,0, param->GetNRow(0)+param->GetNRow(36) );
 252     //
 253     hisQ3D[ith]->SetDirectory(0);
 254     hisQ1D[ith]->SetDirectory(0);
 255     hisQ1D[ith]->SetDirectory(0);
 256     hisQ3DROC[ith]->SetDirectory(0);
 257     hisQ1DROC[ith]->SetDirectory(0);
 258     hisQ1DROC[ith]->SetDirectory(0);
 259   }
 260   //
 261   //
 262   AliRawReader *reader = new AliRawReaderRoot(fileName);
 263   reader->Reset();
 264   AliAltroRawStream* stream = new AliAltroRawStream(reader);
 265   stream->SelectRawData("TPC");
 266   Int_t evtnr=0;
 267   Int_t chunkNr=0;
 268   //
 269
 270   while (reader->NextEvent()) {
 271     Double_t shiftZ= gRandom->Rndm()*250;
 272     //
 273     if(evtnr>=maxEvents) {
 274       chunkNr++;
 275       evtnr=0;
 276       pcstream->GetFile()->mkdir(Form("Chunk%d",chunkNr));
 277       pcstream->GetFile()->cd(Form("Chunk%d",chunkNr));
 278       for (Int_t ith=0; ith<3; ith++){
 279         hisQ1D[ith]->Write(Form("His1DDrift_%d",ith));
 280         hisQ2D[ith]->Write(Form("His2DDrift_%d",ith));
 281         hisQ3D[ith]->Write(Form("His3DDrift_%d",ith));
 282         hisQ1DROC[ith]->Write(Form("His1DROC_%d",ith));
 283         hisQ2DROC[ith]->Write(Form("His2DROC_%d",ith));
 284         hisQ3DROC[ith]->Write(Form("His3DROC_%d",ith));
 285         if (ith<3) continue;
 286         hisQ1D[ith]->Reset();
 287         hisQ2D[ith]->Reset();
 288         hisQ3D[ith]->Reset();
 289         hisQ1DROC[ith]->Reset();
 290         hisQ2DROC[ith]->Reset();
 291         hisQ3DROC[ith]->Reset();
 292       }
 293     }
 294     cout<<"Chunk=\t"<<chunkNr<<"\tEvt=\t"<<evtnr<<endl;
 295     evtnr++;
 296     AliSysInfo::AddStamp(Form("Event%d",evtnr),evtnr);
 297     AliTPCRawStreamV3 input(reader,(AliAltroMapping**)mapping);
 298     //
 299     while (input.NextDDL()){
 300       Int_t sector = input.GetSector();
 301       while ( input.NextChannel() ) {
 302         Int_t    row    = input.GetRow();
 303         Int_t rowAbs    = row+ ((sector<36) ? 0: param->GetNRow(0));
 304         Int_t    pad    = input.GetPad();
 305         Int_t    nPads   = param->GetNPads(sector,row);
 306         Double_t localX  = param->GetPadRowRadii(sector,row);
 307         Double_t localY  = (pad-nPads)*param->GetPadPitchWidth(sector);
 308         Double_t localPhi= TMath::ATan2(localY,localX);
 309         Double_t phi     = TMath::Pi()*((sector%18)+0.5)/9+localPhi;
 310         if (sector%36>=18) phi+=TMath::TwoPi();
 311         Double_t padLength=param->GetPadPitchLength(sector,row);
 312         //  Double_t padWidth=param->GetPadPitchWidth(sector);
 313         //  Double_t zLength=param->GetZLength();
 314         Double_t deltaPhi=hisQ3D[0]->GetXaxis()->GetBinWidth(0);
 315         Double_t deltaZ= hisQ3D[0]->GetZaxis()->GetBinWidth(0);
 316         //
 317         while ( input.NextBunch() ){
 318           Int_t  startTbin    = (Int_t)input.GetStartTimeBin();
 319           Int_t  bunchlength  = (Int_t)input.GetBunchLength();
 320           const UShort_t *sig = input.GetSignals();
 321           Int_t aboveTh[3]={0};
 322           for (Int_t i=0; i<bunchlength; i++){
 323             for (Int_t ith=0; ith<3; ith++){
 324               if (sig[i]>(ith*2)+2) aboveTh[ith]++;
 325             }
 326           }
 327           for (Int_t ith=0; ith<3; ith++){
 328             if (aboveTh[ith%3]>1){
 329               for (Int_t i=0; i<bunchlength; i++){
 330                 //
 331                 // normalization
 332                 //
 333                 Double_t volume3D  = padLength*(localX*deltaPhi)*deltaZ;
 334                 Double_t zIonDrift = (param->GetZLength()-startTbin*param->GetZWidth());
 335                 zIonDrift+=shiftZ;
 336                 if (TMath::Abs(zIonDrift)<param->GetZLength()){
 337                   if ((sector%36)>=18) zIonDrift*=-1;   // c side has opposite sign
 338                   if (sector%36<18) hisQ1D[ith]->Fill(rowAbs, sig[i]/volume3D);
 339                   hisQ2D[ith]->Fill(phi,rowAbs, sig[i]/volume3D);
 340                   hisQ3D[ith]->Fill(phi,rowAbs,zIonDrift,sig[i]/volume3D);
 341                 }
 342                 //
 343                 Double_t zIonROC = ((sector%36)<18)? shiftZ: -shiftZ;  // z position of the "ion disc" -  A side C side opposite sign
 344                 if (sector%36<18) hisQ1DROC[ith]->Fill(rowAbs, sig[i]/volume3D);
 345                 hisQ2DROC[ith]->Fill(phi,rowAbs, sig[i]/volume3D);
 346                 hisQ3DROC[ith]->Fill(phi,rowAbs,zIonROC,sig[i]/volume3D);
 347               }
 348             }
 349           }
 350         }
 351       }
 352     }
 353   }
 354   timer.Print();
 355   delete pcstream;
 356   return 0;
 357 }