PMD/AliPMDHitsRead.C

   1 //
   2 // This macro reads the Hits Tree
   3 //
   4 void AliPMDHitsRead(Int_t nevt = 1)
   5 {
   6
   7   TStopwatch timer;
   8   timer.Start();
   9   TH2F *h2 = new TH2F("h2"," Y vs. X",200,-100.,100.,200,-100.,100.);
  10   //  FILE *fpw = fopen("alipmdhits.dat","w");
  11
  12   AliRunLoader *fRunLoader = AliRunLoader::Open("galice.root");
  13
  14   if (!fRunLoader)
  15    {
  16      printf("Can not open session for file ");
  17    }
  18
  19   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
  20   if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
  21   if (!fRunLoader->TreeK()) fRunLoader->LoadKinematics();
  22
  23   gAlice = fRunLoader->GetAliRun();
  24
  25   if (gAlice)
  26     {
  27       printf("Alirun object found\n");
  28     }
  29   else
  30     {
  31       printf("Could not found Alirun object\n");
  32     }
  33
  34   fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
  35
  36   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
  37
  38   if (fPMDLoader == 0x0)
  39     {
  40       printf("Can not find PMDLoader\n");
  41     }
  42
  43
  44   fPMDLoader->LoadHits("READ");
  45
  46   // This reads the PMD Hits tree and assigns the right track number
  47   // to a cell and stores in the summable digits tree
  48   //
  49
  50   const Int_t kPi0 = 111;
  51   const Int_t kGamma = 22;
  52   Int_t   npmd;
  53   Int_t   trackno;
  54   Int_t   smnumber;
  55   Int_t   trackpid;
  56   Int_t   mtrackno;
  57   Int_t   mtrackpid;
  58   Int_t   xpad = -1, ypad = -1;
  59   Float_t edep;
  60   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
  61   Float_t xPos, yPos, zPos;
  62   Float_t xx, yy;
  63
  64   AliPMDUtility cc;
  65
  66   for (Int_t ievt = 0; ievt < nevt; ievt++)
  67     {
  68
  69       printf("Event Number = %d\n",ievt);
  70       Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
  71       printf("Number of Particles = %d\n",nparticles);
  72       fRunLoader->GetEvent(ievt);
  73       // ------------------------------------------------------- //
  74       // Pointer to specific detector hits.
  75       // Get pointers to Alice detectors and Hits containers
  76
  77       TTree* treeH = fPMDLoader->TreeH();
  78
  79       Int_t ntracks    = (Int_t) treeH->GetEntries();
  80       printf("Number of Tracks in the TreeH = %d\n", ntracks);
  81
  82
  83       TClonesArray* hits = 0;
  84       if (fPMD) hits = fPMD->Hits();
  85
  86       // Start loop on tracks in the hits containers
  87
  88       for (Int_t track=0; track<ntracks;track++)
  89         {
  90           gAlice->GetMCApp()->ResetHits();
  91           treeH->GetEvent(track);
  92           if (fPMD)
  93             {
  94               npmd = hits->GetEntriesFast();
  95               for (int ipmd = 0; ipmd < npmd; ipmd++)
  96                 {
  97                   fPMDHit = (AliPMDhit*) hits->UncheckedAt(ipmd);
  98                   trackno = fPMDHit->GetTrack();
  99
 100                   //fprintf(fpw,"trackno = %d\n",trackno);
 101
 102                   //  get kinematics of the particles
 103
 104                   TParticle* mparticle = gAlice->GetMCApp()->Particle(trackno);
 105                   trackpid  = mparticle->GetPdgCode();
 106
 107                   Int_t igatr = -999;
 108                   Int_t ichtr = -999;
 109                   Int_t igapid = -999;
 110                   Int_t imo;
 111                   Int_t igen = 0;
 112                   Int_t idmo = -999;
 113
 114                   Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
 115                   if (mparticle->GetFirstMother() == -1)
 116                     {
 117                       tracknoOld  = trackno;
 118                       trackpidOld = trackpid;
 119                       statusOld   = -1;
 120
 121                       vx = mparticle->Vx();
 122                       vy = mparticle->Vy();
 123                       vz = mparticle->Vz();
 124
 125                       //fprintf(fpw,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, -1, trackpid, vx, vy, vz);
 126
 127                     }
 128                   Int_t igstatus = 0;
 129                   while((imo = mparticle->GetFirstMother()) >= 0)
 130                     {
 131                       igen++;
 132
 133                       mparticle =  gAlice->GetMCApp()->Particle(imo);
 134                       idmo = mparticle->GetPdgCode();
 135
 136                       vx = mparticle->Vx();
 137                       vy = mparticle->Vy();
 138                       vz = mparticle->Vz();
 139
 140                       //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
 141                       //fprintf(fpw,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
 142
 143                       if ((idmo == kGamma || idmo == -11 || idmo == 11) && vx == 0. && vy == 0. && vz == 0.)
 144                         {
 145                           igatr = imo;
 146                           igapid = idmo;
 147                           igstatus = 1;
 148                         }
 149                       if(igstatus == 0)
 150                         {
 151                           if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
 152                             {
 153                               igatr = imo;
 154                               igapid = idmo;
 155                             }
 156                         }
 157                       ichtr = imo;
 158                     } // end of while loop
 159
 160                   if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
 161                     {
 162                       mtrackno = igatr;
 163                       mtrackpid = igapid;
 164                     }
 165                   else
 166                     {
 167                       mtrackno  = ichtr;
 168                       mtrackpid = idmo;
 169                     }
 170                   if (statusOld == -1)
 171                     {
 172                       mtrackno  = tracknoOld;
 173                       mtrackpid = trackpidOld;
 174                     }
 175
 176                   //printf("mtrackno = %d  mtrackpid = %d\n",mtrackno,mtrackpid);
 177
 178                   xPos = fPMDHit->X();
 179                   yPos = fPMDHit->Y();
 180                   zPos = fPMDHit->Z();
 181
 182                   Float_t time = fPMDHit->GetTime();
 183
 184                   printf("++++++++++ time = %f\n",time);
 185
 186                   edep       = fPMDHit->GetEnergy();
 187                   Int_t vol1 = fPMDHit->GetVolume(1); // Column
 188                   Int_t vol2 = fPMDHit->GetVolume(2); // Row
 189
 190                   Int_t vol3 = fPMDHit->GetVolume(4); // UnitModule
 191
 192
 193                   // -----------------------------------------//
 194                   // For Super Module 1 & 2                   //
 195                   //  nrow = 96, ncol = 48                    //
 196                   // For Super Module 3 & 4                   //
 197                   //  nrow = 48, ncol = 96                    //
 198                   // -----------------------------------------//
 199
 200                   if (vol3 < 24)
 201                     {
 202                       smnumber = vol3;
 203                     }
 204                   else
 205                     {
 206                       smnumber = vol3 - 24;
 207                     }
 208
 209                   xpad = vol1;
 210                   ypad = vol2;
 211
 212                   if(zPos > 361.5)
 213                     {
 214                       cc.RectGeomCellPos(smnumber,xpad,ypad,xx,yy);
 215                       h2->Fill(xx,yy);
 216                     }
 217
 218                 }
 219             }
 220         } // Track Loop ended
 221
 222     }
 223
 224   h2->Draw();
 225
 226   fRunLoader->UnloadgAlice();
 227   fRunLoader->UnloadHeader();
 228   fRunLoader->UnloadKinematics();
 229   fPMDLoader->UnloadHits();
 230
 231   timer.Stop();
 232   timer.Print();
 233
 234 }